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Finite Size Scaling analysis of the avalanches in the 3d Gaussian Random Field Ising 

Model with metastable dynamics 
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Dept. d'Estructura i Constituents de la Materia, Universitat de Barcelona 
Diagonal 647, Facultat de Fisica, 08028 Barcelona, Catalonia 
(Dated: February 1, 2008) 

A numerical study is presented of the 3d Gaussian Random Field Ising Model at T — driven by an 
external field. Standard synchronous relaxation dynamics is employed to obtain the magnetization 
versus field hysteresis loops. The focus is on the analysis of the number and size distribution of 
the magnetization avalanches. They are classified as being non-spanning, ld-spanning, 2d-spanning 
or 3d-spanning depending on whether or not they span the whole lattice in the different space 
directions. Moreover, finite-size scaling analysis enables identification of two different types of non- 
spanning avalanches (critical and non-critical) and two different types of 3d-spanning avalanches 
(critical and subcritical) , whose numbers increase with L as a power-law with different exponents. 
We conclude by giving a scenario for avalanche behaviour in the thermodynamic limit. 

PACS numbers: 75.60.Ej, 05.70.Jk, 75.40.Mg, 75.50.Lk 



I. INTRODUCTION 

Systems with first-order phase transitions exhibit a 
discontinuous change of their properties when driven 
through the transition point. Sometimes, due to the ex- 
istence of energy barriers larger than thermal fluctua- 
tions, such systems evolve following a path of metastable 
states and exhibit hysteresis. Metastable phenomena de- 
velop more often in the case of systems at low temper- 
ature and with quenched disorder. In many cases the 
first-order phase transition occurs, instead of at a cer- 
tain transition point, in a broad range of the driving pa- 
rameter and the discontinuity is split into a sequence of 
jumps or avalanches between metastable states. More- 
over, under certain conditions such avalanches do not 
show any characteristic spatial or time scale: the distri- 
bution of their size and duration becomes a power-law. 
This framework, which has sometimes been called fluc- 
tuationless first-order phase transitions p], Q , is one of 
the basic mechanisms responsible for power-laws in na- 
ture [3 • Experimental examples have been found in a 
broad set ofphysical systems: magnetic transitions 0, 
adsorption gl, superconductivity |6j, martensitic trans- 
formations pi, etc. 

A paradigmatic model for such fluctuationless first- 
order phase transitions in disordered system is the Gaus- 
sian Random Field Ising Model (GRFIM) at T = 
driven by an external field H. The amount of quenched 
disorder is controlled by the standard deviation a of the 
Gaussian distribution of independent random fields act- 
ing on each spin. Metastable evolution is obtained with 
appropriate local relaxation dynamics which assumes a 
separation of time scales between the driving field rate 
dH/ dt and the avalanche duration. The response of the 
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system to the driving field can be followed by measuring 
the total magnetization m(H). The response exhibits the 
above-mentioned metastable phenomena: hysteresis and 
avalanches. 

Since the model was introduced some years ago 0, |j| , 
different studies (numerical and analytical) have been 
carried out in order to characterize the hysteresis loops 
m(H) and the magnetization avalanches jloL HTl Il2l 
fLU ITU . Two of the most well-studied properties are the 
number of avalanches N(cr) and the distribution D(s; a) 
of avalanche sizes s along half a hysteresis loop. For large 
amounts of disorder [a > a c ) the loops look smooth and 
continuous. They consist of a sequence of a large number 
of tiny avalanches whose size distribution D(s; a > a c ) 
decays exponentially with s. On the other hand, for small 
amounts of disorder (p < <r c ), besides a certain number 
of small avalanches, one or several large avalanches pro- 
duce a discontinuity Am in the hysteresis loop. For an 
intermediate critical value a c the distribution D(s, a c ) of 
avalanche sizes s can be approximated by a power-law: 
D{s;a c ) ~s-\ 

Many of the properties of the GRFIM have been un- 
derstood by assuming the existence of a T = critical 
point (a c ,H c ) on the metastable phase diagram. The 
more recent estimation renders: a c = 2.16 ± 0.03 
and H c — 1.435 ± 0.004. Although partial agreement 
on the values of the critical exponents has been reached, 
other features are still controversial. 

One of the fundamental problems is the definition of 
the order parameter. From a thermodynamic point of 
view the discontinuity of the hysteresis loop Am seems 
to be an appropriate order parameter if Am > for 
(7 < <7 C and Am = for a > a c . Nevertheless, in the 
T = numerical simulations, due to the finite size of 
the system and for a given realization of disorder, all the 
magnetization changes are discontinuous. Note that this 
does not occur for standard thermal numerical simula- 
tions in which, due to thermal averaging, magnetization 
is continuous for finite systems. Only finite-size scaling 
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analysis will reveal which are the "large avalanches" and 
whether or not avalanches become vanishingly small in 
the thermodynamic limit. It is thus very important to 
study the properties of the "spanning" avalanches. These 
are avalanches that, for a finite system with periodic 
boundary conditions, cross the system from one side to 
another. In particular it would be interesting to measure 
the number N s (a) of spanning avalanches and their size 
distribution D s (s; a). 

A second unsolved question, related to the previous 
one, is the spatial structure of the avalanches. It has 
been suggested that they are not compact 0, 0]. A 
fractal dimension (df = 1/0.34 < 3) has been estimated 
from the avalanche size distribution [h]. It would be 
interesting to understand how such a fractal behaviour 
may, in the thermodynamic limit, represent a magneti- 
zation discontinuity. 

A third problem is the definition of the scaling vari- 
ables in order to characterize the critical properties close 
to the critical point (a c ,H c ). When focusing on the 
study of avalanche properties, it should be pointed out 
that the scaling analysis is performed by using quan- 
tities (N(a) and D(s,a)) measured recording all the 
avalanches along half a hysteresis loop. The measure- 
ment of non-integrated distributions, i.e. around a cer- 
tain value of H, will require large amounts computing 
effort in order to reach good statistics for large enough 
systems. Therefore, the dependence on the field H is in- 
tegrated out and the distance to the critical point a c is 
measured by a single scaling variable u(a). Although in 
pioneering papers @, El the most usual scaling variable 
u\ = (a — cr c ) /a c was used in order to scale the avalanche 
size distribution, forthcoming studies 0,0, El changed 
the definition to 113 = (a — u c ) /a. Apparently both def- 
initions are equivalently close to the critical point, but 
it can be checked that the "phenomenological" scaling of 
the distributions D(s; a) using w 3 (with u 3 > 0.04) as 
suggested in the inset of Fig. 1 in Ref. 10, is not possible 
when using u\. 

Finite-size scaling analysis has been carried out 
for the number of spanning avalanches N s (a\L). Never- 
theless, such finite-size scaling has not been presented ei- 
ther for the avalanche size distributions D(s;a,L) or for 
the number of non-spanning avalanches N ns (a; L). Most 
of the studies 0, have proposed collapses by neglect- 
ing the fact that simulated systems are finite. There is an 
exception 0, 0] for which the scaling of the avalanche 
distributions with L has been studied. In this case, nev- 
ertheless, the dependence on the distance to the critical 
point has been neglected and, consequently, parameter- 
dependent exponents have been obtained. In our opinion, 
scaling of the avalanche distribution must be studied on a 
two dimensional plane, including a scaling variable which 
accounts for the finite-size L and another which accounts 
for the distance to the critical point. 

Previous studies have provided simulations of very 
large system sizes (up to L = 1000) ^4|. This has been 
advantageous for the study of self-averaging quantities. 



Nevertheless, the properties of the spanning avalanches 
are non-self averaging. This is because, as will be shown, 
the number of spanning avalanches per loop does not 
grow as L 3 . This means that, in order to obtain better 
accuracy, it is more important to perform averages over 
different disorder configurations (which will be indicated 
by (•) ) than to simulate very large system sizes. 

In this paper we present intensive numerical studies of 
the metastable 3d-GRFIM and focus on analysis of the 
spanning avalanches. In section [n] the model, the defini- 
tion of a spanning avalanche and the details of the numer- 
ical simulations are presented. In section ITTT1 raw numeri- 
cal results are given. In section llVl some of the Renormal- 
ization Group (RG) ideas will be reviewed, which will be 
taken into account for the analysis of the critical point. 
A finite-size scaling analysis of the avalanche numbers is 
presented in section [V] The same analysis for size dis- 
tributions and their fc-moments are presented in sections 
IVII and IVlII respectively. Section IVIIII presents a discus- 
sion on the behaviour of magnetization. The discussion 
of the results in relation with previous works is presented 
in section IIXI Finally in section ^ a full summary and 
conclusions is given. 

II. MODEL 

The 3d-GRFIM is defined on a cubic lattice of size 
L x L x L. On each lattice site (i — 1, . . . , L 3 ) there is a 
spin variable Si taking values ±1. The Hamiltonian is: 

71. n. L 3 L 3 

H = -J2S i Sj-J2h i S i -H^2S i , (1) 

i,j i—1 i—1 

where the first sum extends over all nearest-neighbour 
(n.n) pairs, H is the external applied field and hi are 
quenched random fields, which are independent and are 
distributed according to a Gaussian probability density 

1 

dP{h l ) = ^=^e~^dh l , (2) 
V2na 

where the standard deviation a, is the parameter that 
controls the amount of disorder in the system. Note that 
(hi) = and (hf) = a 2 . 

The system is driven at T — by the external field 
H. For H = +00 the state of the system which mini- 
mizes Ti. is the state with maximum magnetization m — 

Ylf=i Si/L 3 = 1. When the external field H is decreased, 
the system evolves following local relaxation dynamics. 
The spins flip according to the sign of the local field: 

6 

hi + H + J2 S 3 > ( 3 ) 
3=1 

where the sum extends over the 6 nearest-neighbouring 
spins of Si. Avalanches occur when a spin flip changes 
the sign of the local field of some of the neighbours. This 
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averaged number 


notation 


avalanches 
spanning avalanches 
non-spanning avalanches 
critical non-spanning avalanches 
non-critical non-spanning avalanches 
ld-spanning avalanches 
2d-spanning avalanches 
3d-spanning avalanches 
critical 3d-spanning avalanches 
subcritical 3d-spanning avalanches 


N(a,L) 
N s (a,L) 
N ns (a,L) 
N nsc (a, L) 
N ns o(&, L) 
Ni (cr, L) 
N 2 (a,L) 
Nz{a, L) 
N 3c (a,L) 
N 3 ~(a,L) 


normalized size distribution 


notation 


avalanches 
spanning avalanches 
non-spanning avalanches 
critical non-spanning avalanches 
non-critical non-spanning avalanches 
ld-spanning avalanches 
2d-spanning avalanches 
3d-spanning avalanches 
critical 3d-spanning avalanches 
subcritical 3d-spanning avalanches 


D{s;<7, L) 
D 3 (s;a,L) 
D ns (s; a, L) 
D nsc (s;a,L) 
D nB o(s; cr, L) 
D 1 (s;a,L) 
D 2 {s;a,L) 
D 3 (s;a,L) 
D 3c (s;a,L) 
D 3 -(s;a,L) 



TABLE I: Notation of the studied quantities in this work. All 
the quantities refer to the analysis of half a hysteresis loop and 
are obtained after averaging over many different realizations 
of disorder. 

may start a sequence of spin flips which occur at a fixed 
value of the external field H , until a new stable situation 
is reached. H is then decreased again. This "adiabatic" 
evolution corresponds to the limit for which avalanches 
are much faster than the decreasing field rate. Note that, 
once the local random fields are fixed, the metastable evo- 
lution is completely deterministic, no inverse avalanches 
may occur and the hysteresis loops exhibit the return 
point memory property 0. 

The size of the avalanche s corresponds to the number 
of spins flipped until a new stable situation is reached. 
Note that the corresponding magnetization change is 
Am = 2s/L 3 . 

For a certain realization of the random fields, corre- 
sponding to a given value of cr, we have recorded the 
sequence of avalanche sizes during half a hysteresis loop, 
i.e. decreasing H from +00 to —00. The two main quan- 
tities (see TableHJ that are measured after averaging over 
different realizations of disorder, are the total number 
of avalanches per loop N(a, L) and the distribution of 
avalanche sizes D(s; cr, L), normalized so that: 

L 3 

^ J D( S ; C r,L) = l . (4) 

s=l 

Note that given this normalization condition and the fact 
that s is a natural number, then D(s; cr, L) < 1. 

The numerical algorithm we have used is the so-called 
brute force algorithm propagating one avalanche at a 
time [l it . We have studied system sizes ranging from 
L = 5(L 3 = 125) to L — 48(L 3 = 110592). The mea- 



Closure relations 


AT AT I AT 

N ns = N nsc + N ns o 

N s = Ni+N 2 + N 3 
N 3 = N 3c + iV 3 - 


Normalization condition 


Es-iA,(s;cr,L) = 1 


Distribution relations 


ND = N s D a + N ns D ns 

N S D S = N1D1 + N 2 D 2 + N 3 D 3 
N 3 D 3 = N 3c D 3c + N 3 -D 3 - 



TABLE II: Main mathematical relationships between the 
quantities defined in Table The dependence on cr, L and s 
has been suppressed in order to clarify the Table. The sub- 
script a stands for all the possible sub-indices in Table Q 



sured properties are always averaged over a large number 
of realizations of the random field configuration for each 
value of cr. Typical averages are performed over a num- 
ber of configurations that ranges between 10 5 for L < 16 
to 2000 for L = 48. 

We have used periodic boundary conditions: the nu- 
merical simulations correspond, in fact, to a periodic in- 
finite system. Therefore, strictly speaking, all avalanches 
are infinite. Nevertheless, we need to identify which 
avalanches will become important in the thermodynamic 
limit. The definition that best matches this idea is the 
concept of spanning avalanches: those avalanches that, 
at least in one of the x, y or z directions, extend over 
the length L. This definition is very easy to imple- 
ment numerically in the brute force algorithm. Span- 
ning avalanches are detected by using three (cc, y, z) mask 
vectors of size L whose elements are set to at the be- 
ginning of each avalanche. During the evolution of the 
avalanche the mask vectors record the shade of the flip- 
ping spins along the three perpendicular directions (by 
changing the 0's to l's). When the avalanche finishes, 
it can be classified as being non-spanning, ld-spanning 
, 2d-spanning or 3rf-spanning depending on the number 
of such mask vectors that have been totally converted 
to 1. The number and size distribution of Id, 2d and 
3d-spanning avalanches is also studied and averaged over 
different realizations of disorder. Table [I] shows the defi- 
nitions of avalanche numbers and distributions that will 
be used throughout the paper. In Table [D] the list of 
mathematical relations between the avalanche numbers 
and distributions is given. We will use the subscript a to 
indicate any of the avalanche numbers or distributions in 
Table 13 

It should be mentioned that, although the definition 
of spanning avalanches used in this paper is equivalent 
to the definition in previous works [l3t fLU Hfil| . the aver- 
age number of spanning avalanches N s , in some cases, 
does not coincide with the previous estimations. We 
guess that the reason is because, in previous works, the 
method used to count spanning avalanches was aver- 
aging twice the 2d-spanning avalanches and was aver- 
aging three times the 3d-spanning avalanches. There- 
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fore, in order to compare, for instance with Ref. \l$L one 
should take into account that their number of spanning 
avalanches N is not equal to the present N s but satisfies: 
N = (Ni + 2N 2 + 3iV 3 )/3. Moreover, we should point 
out the following remark before presenting the data. As 
a consequence of the numerical analysis, several "kinds" 
of avalanches will be identified (see Table P). Such a 
separation in different kinds will, in some cases, be justi- 
fied by the measurement of different physical properties 
(such as whether the avalanche spans the lattice or not) 
but, in other cases, will be an "a priori" phenomcnolog- 
ical hypothesis to reach a good description of the data. 
Although some authors will prefer to identify such new 
"kinds" of avalanches as "corrections to scaling", it will 
turn out that after the finite size scaling analysis we will 
be able to identify which different physical properties 
characterize each "kind" of avalanche. 



III. NUMERICAL RESULTS 

Fig. n shows an example of the distribution of 
avalanches D(s; a, L) on a log-log scale for three values of 
a corresponding to a system with size L = 24. The qual- 
itative behaviour of D(s;a,L) is that already described 
in the introduction: when a is decreased the distribution 
changes from being approximately exponentially damped 
(a > er c ) to a distribution exhibiting a peak for large val- 
ues of s (cr < <t c ). Therefore, one can suggest that at 
the critical value a c the distribution exhibits power-law 
behaviour. Nevertheless, it is also evident from Fig. 2] 
that the finite size of the system masks this excessively 
simplistic description. Only after convenient finite-size 
scaling analysis shall we discover which features remain 
in the thermodynamic limit. 

The peak occurring for a < a c is basically caused by 
the existence of spanning avalanches. This is shown in 
Fig.Hwhere the peak in D(s; a, L) [Fig.^a)] is compared 
with the two contributions D ns {s;a,L) and D s (s; a, L) 

pig. EDO] • 

As can be seen the distribution of spanning avalanches 
D s (s; a, L) is far from simple. It exhibits a multipeak 
structure caused by the contributions from Di(s;a,L) , 
D2(s;a,L) and D 3 (s;<r,L) shown in Fig. 0;. Moreover, 
Dz{s\cr, L) itself also exhibits two peaks suggesting that 
the 3d-spanning avalanches may be of two different kinds. 
We shall denote critical 3d-spanning avalanches (indi- 
cated by the subscript 3c) as those corresponding to the 
peak on the left and subcritical 3d-spanning avalanches 
(indicated by the subscript 3—) as those corresponding to 
the peak on the right. As will be explained below, the Id- 
spanning avalanches, the 2d-spanning avalanches and the 
critical 3d-spanning avalanches do not exist in the ther- 
modynamic limit except when a = a c . This is the reason 
for having chosen the word "critical" for this kind of 3d 
spanning avalanche. It will also be shown that, in the 
thermodynamic limit, subcritical 3d-spanning avalanches 
only exist for a < a c . As regards the non-spanning 
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FIG. 1: Avalanche size distribution corresponding to (a) 
a — 1.7, (b) a = 2.21 and (c) a — 2.6. Data have been 
obtained from a system with size L — 24 after averaging over 
10 J realizations of the disorder. 



avalanches, they will also be classified into two types at 
the end of this section, although this separation cannot 
be deduced from the behaviour in Fig.J^b). 

Fig. |21 shows the evolution of D\(s, a, L), D 2 {s. a, L) 
and Dz(s 1 a, L) when a is increased. Note that the right- 
hand peak of D3(s;a,L) shifts to smaller values of s 
and becomes flat, indicating that the mean size of these 
subcritical 3d-spanning avalanches decreases. Moreover, 
above a c [Fig.|3Jd)] this right-hand peak disappears and 
a peak on the left emerges. 

Besides the normalized distributions, it is also inter- 
esting to analyze the actual average numbers of spanning 
avalanches Ni(<j,L), N2{u,L) and Ns(a,L), which also 
exhibit singular behaviour at a c as shown in Fig. 0] 

From the direct extrapolation of the data correspond- 
ing to different system sizes to L — > oo, we can make 
the following assumptions: in the thermodynamic limit 
Ni(a) and A^er) will display a S- function discontinuity 
at cr c . Nsla) will display step-like behaviour: for a < a c 
there is only one 3d-spanning avalanche, for a > a c 
there are no 3d-spanning avalanches and at a = a c the 
data supports the assumption that N3 will also display 
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FIG. 2: Analysis of the different contributions to D(s;a,L) 
for a — 2.317 and L = 16. Data corresponds to an average of 
2 x 10 realizations, (a) Full distribution, (b) distribution of 
spanning avalanches D a (s; a, L); (c) Distributions Di(s; a, L), 
D2{s;a,L) and Dz{s;o, L). 



a ^-function singularity at the edge of the step function. 
This reinforces the suggestion that there are two differ- 
ent types of 3d-spanning avalanches: as will be shown, 
in the thermodynamic limit, the number of subcritical 
3d spanning avalanches iV 3 _ behaves as a step function, 
whereas the number of critical avalanches N^ c exhibits 
divergence at a c . 

The total number of spanning avalanches N s (cr, L) 
and non-spanning avalanches N ns (a, L), are displayed in 
Figs. EI a) and 02b) respectively. N s (a,L) shows, as a 
result of the divergence of N^, Ni and N2, a (S-function 
singularity at <r c when L — > 00 suggesting that the criti- 
cal point is characterized by the existence of 00 spanning 
avalanches. We would like to point out that previous 
studies have not clarified this result for the 3d-GRFlM 


The analysis of N ns is more intricate. Fig. OJb) shows 
that N ns (a,L) grows with a and L. For large amounts 
of disorder (a — > 00) one expects that the hysteresis loop 
consists of a sequence of non-spanning avalanches of size 
1. Therefore, their number will equal L 3 . To reveal this 
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FIG. 3: Analysis of the dependence of Di(s,a; L) (top), 
D2(s,a; L) (middle) and Ds(s, a; L) (bottom) with a. Data 
correspond to averages of 2 x 10 5 realizations of a L — 10 
system with (a) a — 1.9, (b) a = 2.2, (c) a = 2.5 and (d) 
a = 2.8. 



behaviour Fig. shows the dependence of N ns (a,L)/L 3 
as a function of a. One expects that these lines tend to 
1 when a — ► 00. Moreover, a closer look reveals that at 
er c ~ 2.21, there is a contribution to N ns (a, L)/L 3 which 
decreases with system size. For low values of a one ex- 
pects that non-spanning avalanches always exist, except 
at a = 0. This last statement can easily be understood by 
noticing that an approximate lower bound to the number 
of non-spanning avalanches can be computed by analyz- 
ing how many of the spins Si will flip by themselves, 
independently of their neighbours, due to the fact that 
the local field hi is either larger than 6 or smaller than 
—6. This analysis renders N ns /L 3 > $ err (6/cr) where 
$ err is the error function. 

From these considerations, we expect that for L — » 00 
the curves in Fig. EJtend to a certain limiting behaviour 
which increases smoothly from to 1. This can also 
be appreciated in the inset in Fig. HJ1 which shows the 
behaviour of N ns (a,L)/L 3 as a function of L for four 
different values of the amount of disorder: a = 1.7, cr = 
2.21 ~ er c , a — 2.5 and a — 2.7. The four curves exhibit 
a tendency to extrapolate to a plateau when L — > 00. 
For the case of a ~ a c an estimation of the extrapolated 
value is N ns (a c , L)/L 3 -> 0.028. 

Consequently, it is necessary to consider the existence 
of, at least, two kinds of non-spanning avalanches. Those 
whose number N ns o increases as L 3 will be denoted as 
non-critical non-spanning avalanches (with the subscript 
nsO), and those whose number N nsc increases with L with 
a smaller exponent will be called critical non-spanning 
avalanches (with the subscript nsc). In fact, a log- log plot 
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FIG. 4: Number of spanning avalanches in Id (a), 2d (b) and 
3d (c) as a function of a. The different curves correspond to 
L = 5, 8, 10, 12, 16, 24, 32 and 48 as indicated by the legend. 
Lines are a guide to the eye. 



of N ns (<r c , L) I L 3 — 0.028 versus L provides an estimation 
for this exponent N nsc (a c ,L) ~ 0.085L 2 02 . 

All the assumptions that have been presented, corre- 
sponding to behaviour in the thermodynamic limit, will 
be confirmed by the finite-size scaling analysis presented 
in the following sections. 



IV. RENORMALIZATION GROUP AND 
SCALING VARIABLES 

The basic hypothesis for the analysis of the above re- 
sults using RG techniques is the existence of a fixed point 
in the multidimensional space of Hamiltonian parame- 
ters. This fixed point sits on a critical surface which 
extends along all the irrelevant directions. By changing 
the two tuneable parameters a and H, the critical sur- 
face can be crossed at the critical point (<7 C , H c ). As has 
been explained in the Introduction concerning the anal- 
ysis of the avalanche number and size distributions, the 
dependence along the external field direction H has been 
integrated out. One expects that such integration may 
distort some of the exponents and the shape of scaling 
functions, but not the possibility of an RG analysis. This 




I i i i i i i i i i i i I 
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FIG. 5: (a) Total number of spanning avalanches N s (a,L) 
and (b) non-spanning avalanches as a function of a for dif- 
ferent system sizes L as indicated by the legend. Lines are a 
guide to the eye. 



is because the integration range crosses the critical sur- 
face where the divergences occur. 

For a L — > oo system we assume a unique scaling vari- 
able u(a) which measures the distance to a c . The depen- 
dence of u on (t should be smooth, but its proper form 
is unknown [l^. We will discuss three different possibil- 
ities: 

1. The standard choice is to use a dimensionless first 
approximation by expanding u(a) as: 



Nevertheless, in general, the correct scaling vari- 
ables may have a different dependence on a. For in- 
stance, this may be due to the existence of other rel- 
evant parameters, such as the external field, which 
has been integrated out. 

2. A second choice is to extend the expansion of u(a) 
to second order by including a fitting amplitude A: 
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FIG. 6: Number of non-spanning avalanches N ns (a,L) di- 
vided by L s as a function of a and different system sizes, as 
indicated by the legend. The inset shows the behaviour of 
the same quantity as a function of L for different values of 
a. The dashed line indicates the value N/L 3 = 0.028 and the 
continuous line is a fit of the behaviour proposed in equation 



FIG. 7: Comparison of the behaviour of the three choices for 
the scaling variable u, discussed in the text. We have used 
o c = 2.21. 



let us consider a renormalization step of a factor b close 
to the critical point 0,113]! such that lengths behave as: 



u 2 = a —^ + A(°—^-) 2 . (6) 

3. A third choice, which has been used in previous 
analyses and may be "phenomenologically" justi- 
fied is: 



Note that the Taylor expansion of this function is: 




L b = b~ l L . (9) 

(The variables with the b subscript correspond to the 
renormalized system). We expect that after re-scaling 
the variable u, measuring the distance between a and a c 
changes as: 

u b = b 1/u u , (10) 

which is the standard definition of the exponent v which 
characterizes the divergence of the correlation length 
when a — ► a c . Under the same renormalization step we 
assume that: 

s b = b- d °s . (11) 



Fig. shows the behaviour of the three scaling vari- 
ables U\ (<r) , u 2 (<r) and U3 (a) . For the representation of 
U2 we have chosen A = —0.2, which is the result that 
we will fit in the following sections. The three choices 
are equivalently close enough to the critical point. Nev- 
ertheless, the amplitude of the critical zone, where the 
scaling relations are valid, may be quite different. Since 
A < 0, the variable ui cannot be used for cr > cr c since 
ui{a) shows a maximum at a — 7.735 = 3.5<7 C . A similar 
problem occurs with U3 since, due to its asymptotic be- 
haviour (1L3 — > 1 for o — ► 00) , systems with a large value 
of a cannot be distinguished one from another. 

For the finite system, the magnitudes presented in Ta- 
ble [I] depend on cr, L and, in the case of the size distri- 
butions, on s. In order to identify the scaling variables, 



This latter equation introduces an exponent d a (which 
has been called 1 /Vex by other authors Q ) and can be in- 
terpreted as the fractal dimension of the avalanches. As 
mentioned in the previous section, we expect to find dif- 
ferent types of avalanches. As will be shown numerically 
from the scaling plots in the following sections, it is pos- 
sible to assume that the different types of avalanches be- 
have with the same fractal dimension d a = df, except for 
subcritical 3d-spanning avalanches (for which <i 3 _ 7^ dt) 
and non-critical non-spanning avalanches. 

Close to the critical point the system exhibits invari- 
ance under re-scaling. Therefore, in order to propose a 
scaling hypothesis of the numbers of avalanches N a and 
the avalanche size distributions D a , it is important to 
construct combinations of the variables u, L and s, which 
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remain invariant after renormalization. We find: 

L b ' u b = L ' u , 

al /vda u b = s 1/vd -u . 



(12) 
(13) 
(14) 



Note that these three invariant quantities are not inde- 
pendent since equation l|12(l corresponds to equation l|14|l 
multiplied by equation l|13|) to the power of —l/vd a . 



V. SCALING OF THE NUMBERS OF 
AVALANCHES N a (a,L) 

The discussion in the previous section, enables us to 
propose the following scaling hypothesis: 



N a {a,L) = L 8a N a (uL 1 /^ 



(15) 



The exponent 9 a characterizes the divergence of the 
avalanche numbers at the critical point when L — > oo. 
Note that this definition of 9 a (which is the same used in 
previous works is not consistent with the standard 
finite-size scaling criterion for which the magnitudes grow 
with exponents divided by v. [T^.l20ll2l| . 

As will be shown, the behaviour of the number of Id- 
spanning avalanches, 2d-spanning avalanches and critical 
3d-spanning avalanches can be described with the same 
value of 8i = 02 = 9 3c = 9, so that: 



N 1 (a,L)=L 8 N 1 (uL l ' v 



(16) 



N 2 {a,L) = L e N 2 (uL 1 ^ , (17) 



N 3c {a,L) = L e N 3c (uL 1 /^ 



(18) 



We have tried, without success, to scale the number of 
critical non-spanning avalanches with the same exponent 
9. We therefore need to define a different exponent 9 nsc , 
so that: 



Nr, 



(19) 



As regards the number of N 3 _ avalanches, which is dif- 
ferent from zero away from the critical point in the ther- 
modynamic limit, we propose a scaling hypothesis that 
is compatible with the limiting behaviour at a — and 
d = oo . This leads us to the following assumption: 



N 3 -{tr,L) = N 3 _ (uL 1 ^ , 



(20) 



since in the absence of disorder we expect that the hys- 
teresis loop displays a single avalanche of size L 3 , and, 



consequently the number of avalanches must be N 3 ^ = 1 
independent of the value of L. 

As regards N ns0 it has already been discussed that such 
avalanches will exist in the thermodynamic limit for all 
values of a. Moreover, they are probably not related to 
critical phenomena at a c . For this reason we propose the 
following non-critical dependence: 



N ns0 {a,L) = L 3 N ns0 (a) 



(21) 



In particular, as already mentioned, for large values of 
disorder (a — > +00) these avalanches will be of size s = 1, 
and their number will be N ns o(oo) = L 3 . 

It should also be mentioned that the scaling equations 
(|15f) admit alternate expressions by extracting the vari- 
able uL~ x l v with the appropriate power so that it cancels 
out the dependence on L: 



N a {a,L) = \u\- ve «N, 



(22) 



Nevertheless, such expressions are not very useful for the 
scaling analysis close to a c since they will display a large 
statistical error due to the fact that u — > when a — * a c . 

Figs. |H1 and El show the best collapses corresponding to 
equations (|16fl and l|17fl with the three different choices 
for the variable u, explained in section II VI Data corre- 
sponding to L = 5, 8, 10, 12, 16, 24, 32 and 48 have been 
used. The quality of the collapses close to a c is quite 
good in the three cases. The values of the free param- 
eters that optimize each collapse are indicated on the 
plots. By visual comparison one can see that u 2 is the 
best choice since it allows the smaller sizes to collapse 
too. Of course, this is because the collapses in this case 
have an extra free-parameter A. As regards the quality 
of the overlaps, no remarkable differences are observed 
between the choices u\ and u 3 . In the following collapses 
we will use u 2 with A = —0.2. Thus, the best estimations 
of the free parameters are: a c — 2.21 ±0.02, v — 1.2 ±0.1 
and 9 = 0.10 ±0.02. 

The procedure for improving the collapse of the data 
corresponding to different system sizes, which will be 
used many times throughout this paper, renders what 
we will call "the best values" of the free parameters. Er- 
ror bars represent the estimated range of values for which 
the collapses are satisfactory. We would like to note that 
the obtained value of a c (for the three choices of the vari- 
able u) is slightly higher than the value a c = 2.16 ± 0.03 
proposed in Ref. 0. 

It is interesting to note that the scaling functions Ni 
and N 2 can be very well approximated with Gaussian 
functions. The fits, shown in Figs. Elb) and|!Jb), have 
three free parameters: the amplitude a, the peak position 
x and the width w. The best numerical estimations are: 
ai = 0.946±0.004, x x = 2.691±0.008, wi = 1.293±0.008, 
a 2 = 0.497 ±0.002, x 2 = 2. 227 ±0.007 and w 2 = 1.086 ± 
0.007. 

From the fact that the scaling functions in Figs. Etb) 
and b) are bounded and go exponentially to zero for 
u 2 L- x l v -> ±00 (as can also be checked from a log-linear 



9 




FIG. 8: (a) Scaling plot of the number of ld-spanning 
avalanches according to equation 116H using ui as the scal- 
ing variable. The values of the free parameters for the best 
collapses are indicated by the legend. Symbols correspond to 
the sizes indicated in the legend of Fig. [I] (b) Same plot, but 
using the scaling variable 112. Note that in this case there is 
an extra free parameter, (c) Same plot but using U3. The 
continuous line in (b) shows the fit of a Gaussian function. 



FIG. 9: (a) Scaling plot of the number of 2d-spanning 
avalanches according to equation 117H using ui as the scal- 
ing variable. The values of the free parameters for the best 
collapses are indicated by the legend. Symbols correspond to 
the sizes indicated in the legend of Fig. [I] (b) Same plot, 
but using the scaling variable U2- Note that in this case there 
is an extra free parameter, (c) Same plot but using U3.The 
continuous line in (b) shows the fit of a Gaussian function. 



plot) one can deduce that, in the thermodynamic limit, 
ld-spanning avalanches and 2d-spanning avalanches only 
exist at a = a c . Their numbers increase as L 010 with 
amplitudes iVi(O) = 0.12 ±0.01 and AT 2 (0) = 0.07±0.01. 
Moreover, the peaks of the scaling functions N\ and N2 
which are displaced from ui = 0, account for the fact 
that for a finite system the maximum number of Id and 
2d spanning avalanches occurs for a certain a c (L) which 
shifts towards a c from above. 

As regards the 3d spanning avalanches, according to 
the previous discussions one must consider the contribu- 
tions from A^3 C and A^_. From the scaling assumptions 
(|18H and l|20|) and the last closure relation in Table ITU one 



can write: 

N 3 (a,L)=L 8 N 3c (uL 1 / u )+N s _(uL 1 ^ . (23) 

This equation indicates that N^a^L) cannot be col- 
lapsed in a straightforward way. We propose here a 
method to separate the two contributions in equation 
(|23|l . This method, which we will call double finite-size 
scaling (DFSS), will be used several times throughout the 
paper for the analysis of similar equations. By choosing 
two systems with sizes Li and L2 and amounts of disor- 
ders u\ and (72 so that u(ai)L\^ u = u^a^)^^ ■> one can 
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write: 



(u{a x )L\' v 



(u{<Ti)L\' v ) 



N 3 _ {u{a 2 )Ll /v ) = (24) 



Ln 



Nsc (u{a 2 )L\ ,u ) = 
N s {ax,L 1 )-N s (a2,L 2 ) 



(25) 



Thus, we can check for the collapse of data corresponding 
to different pairs of (Lx^L^. From the numerical point 
of view, the DFSS method works quite well. An analysis 
of error propagation reveals that the scaling function cor- 
responding to the contribution with a smaller exponent 
will display more statistical errors. 

Fig. EH shows the results of the DFSS analysis of N 3 
according to equation l|23|) . The different symbols, in 
this case, indicate the values of L\ and L2 used for each 
set of data. Fig. HUf a) corresponds to N^-i^L 1 /") and 
Fig. rrUT b) corresponds to N^^L 1 ^). It should be em- 
phasised that such collapses are obtained without any 
free parameter. The values of 9, a c , v and A are taken 
from the previous collapses of N\ and N% ■ 

Again, from the shape of the scaling functions we 
can deduce the behaviour in the thermodynamic limit: 
from the crossing points of the scaling functions with the 
u 2 = axis, we find that A 3c (cr c , L) = (0.16 ± 0.02)L 010 
and N 3 ^(cr c ,L) = 0.79 ± 0.02. As occurred previously 
with the number of 2d and Id spanning avalanches, N 3c 
can also be very well approximated with a Gaussian func- 
tion with amplitude a 3c = 0.706 ± 0.005, peak position 
x 3c = 1.244 ± 0.007 and width w 3c = 0.802 ± 0.009. 
The fact that N 3c {u2L 1 / v ) vanishes exponentially for 
U2L X I V — » ±00 confirms that, in the thermodynamic 
limit, such avalanches only exist at the critical point. 
Furthermore, from the fact that A 3 _ tends to 1 and 
to exponentially fast when U2L 1 / 1 ' — > ±00 we deduce 
that one subcritical 3d-spanning avalanche will exist for 
a < <7 C and there will be none above this value. 

To end with the analysis of the number of avalanches, 
we will separate the two contributions to N ns : 

N ns (a, L) = N nsc (uL 1 '^ + L 3 N ns0 (a) (26) 

In this case the DFSS method cannot be applied since 
N nsc and N ns o depend on different variables. A first 
check of the validity of this hypothesis has already been 
presented in section III. The fit of equation l|26|l to the 
data corresponding to a — a c (u — 0) , shown in the inset 
of Fig. gives estimations of 6 nsc ~ 2.02 , N n sO (c c ) = 
0.028, and N nsc (0) = 0.085. Furthermore, we can also 
check that the derivative with respect to a behaves as: 



1_ dN ns (a, L) 
L 3 



dc 



4-3 



K s0 (<Tc) 

(27) 
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FIG. 10: (a) Numerical estimation of N 3 - {u 2 L 1/,J ) and (b) of 
N3c(u2L 1 i l "). Data have been obtained according to equations 
1241 and 1251 Symbols, according to the legend, indicate the 
values of L\ and L2 used for obtaining each data set. The 
continuous line in (b) is a fit of a Gaussian function. 



Fig- II If a) demonstrates that the data (estimated us- 
ing a two-point derivative formula) is compatible with 
this behaviour. The line shows the best fit (with two 
free parameters: N' nsc (0) and N^sO (°c)) of the func- 
tion 127| l with 9 nsc + \ jv — 3 = —0.15. One obtains 
N' nsc (0) = -0.136±0.011and N' ns0 (cr c ) = 0.102±0.003. 
The good agreement is a test of the dependence with the 
variables uL 1 / 1 ' and a of the functions N nsc and N n sO 
respectively. To go further into the analysis of N ns , one 
must provide some extra hypothesis on the shape of the 
scaling functions. Given the fact that we have found al- 
most a perfect Gaussian dependence of the scaling func- 
tions Ni, N2 and N 3c one can guess that N nsc will also 
have a Gaussian dependence. By forcing the Gaussian 
function to satisfy N nsc (0) — 0.085 and the fact that 
that N^ SC (Q) = —0.136 (from previous estimations) we 
end up with a trial function with a single free parame- 
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ter that should be enough to satisfactorily scale the data 
from Fig. [S] 
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FIG. 11: (a) Analysis of the consistency of Eq. J23- The 
points correspond to numerical data and the line is the best 
fit (with two free parameters) by fixing the value of nac + l/v 
to the previous estimations, (b) Scaling plot revealing the be- 
haviour of Nnso(o')- The continuous line shows the behaviour 
of the approximate lower bound $ err (6/cr) estimated from 
the spins that flip independently of their neighbours. The 
inset shows the Gaussian function N nS c{u2L 1 ^ v ) used for the 
separation of the two contributions to N ns . 

The best collapse is shown in Fig. IllT b) which cor- 
responds to N ns o(cr). The function N nsc used for the 
collapse is shown in the inset and corresponds to a Gaus- 
sian function with amplitude a nsc — 0.085, peak position 
x n sc — —0.6 and width w nsc = 1.485. It is interesting to 
note that the peak position of this scaling function occurs 
at a value U2L 1 / 1 ' = x nsc < as opposed to the case of 
the previous scaling functions N±, N2 and AT 3c for which 
the peak position was at u^L 1 ^ > 0. This indicates that 
the properties of the Id, 2d and 3c critical avalanches 
have opposite shifts with finite size L compared to the 
nsc critical avalanches. 

To end with the analysis of the number of non- 
spanning avalanches it is interesting to compare the 



function N ns o(a) with the approximate lower bound 
err (6 /a)) discussed in section III, which is represented 
by a continuous line in Fig. Illf b). The difference be- 
tween the two curves, which becomes bigger when <j 
increases, is due to the existence of clusters of several 
spins (not considered in the extremely facile analysis pre- 
sented here) that flip independently of their neighbours 
contributing to the number of non-critical non-spanning 
avalanches. 



VI. SCALING OF THE DISTRIBUTIONS OF 
SIZES D a (s;a,L) 

Close to the critical point there are different ways to ex- 
press the invariance of the size distributions correspond- 
ing to different choices of a pair of invariants among the 
three exponents proposed in equations 112(1. (|13fl and 
(|14fl . For any generic distribution D a (s;a, L) one can 
write the following nine generic expressions: 



D a (s-a,L) = L- T " d "b a (sL- d ",uL 1 / u ^ , (28) 

D a (s;a,L) = L- T ^D a {s^ vd ''u,uL 1 l v \ , (29) 

D a (s;a,L) = L Tada D a i^s 1 ^ udc "u, sL~ dc "^j , (30) 

D a (s;a,L) = s~ T ° D a (sL~ d = ,uL 1/v ^ , (31) 

D a (s;a,L) = s^A* (sV vda u,uL 1 / u ^ , (32) 

D a (s;a,L) = s - T ^b a {s 1/vd "u,sL~ d ^ , (33) 

D a (s;a,L) = \u\ T <*'" l <*D cc (si"^, uL 1 ^ , (34) 
D a (s;a,L) = \u\ TaVd "D a (s 1 / ud «u,uL 1 ^ , (35) 

D a (s;a,L) = \u\ T « vd « b a (s 1/vd «u, sL~ d ^ . (36) 

Although we have used the generic index a, it is evident 
that such assumptions can only be proposed for the dis- 
tributions of avalanches of a single kind, i.e. D±, D2, 
Z?3 C , £>3_ and D nsc . For the composite distributions D3, 
D s , D ns and D, one expects mixed behaviour, and con- 
cerning D nS Q we cannot expect a dependence on uL 1 l lJ . 
The exponents r Q could also be different for the differ- 
ent kinds of avalanches, but as will be discussed in the 
following paragraphs, in all cases T a = 1 except for r„ sc , 
which will take a larger value. 

As argued before, when scaling the numbers of 
avalanches, the last three expressions (13411 . (|35|) and 136|) 
are not very useful for the numerical collapses because 
they introduce large statistical errors. Moreover, when 
trying to check the collapses expressed by equations lt2^|) 
and (|32|l . the two independent variables of the scaling 
function converge to zero when the critical point is ap- 
proached. Thus, such a collapse cannot be checked for 
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u = 0. Therefore, the interesting scaling equations are 
(HEJ, EU and 

The behaviour of the scaling functions is, in some cases, 
restricted by the normalization conditions. If scaling 
holds for the whole range of s = 1, • • • , L 3 , from equation 
(|28J) . one can write: 

L 3 

J2 L ' Tada D a (sL-d^uL 1 /^ = 1 . (37) 

s=l 

If < d a < 3, by defining a new variable x = sL~ dc ", 
for large L, the above expression is transformed into the 
following integral: 

7r (r =- 1)d " J dxD a (x, uL 1 '^ = 1 . (38) 

For those distributions for which the integral converges, 
it is necessary that T a = 1. We expect that this condition 
can be applied to the cases of Di, D 2 , D 3c and -D3-. In 
these four cases, as can be seen in Fig- Etc), the distri- 
butions show a marked decay in the two limits of s — > 
and s — > L 3 . (Note that the plots have logarithmic scales 
and that D 3c and D3- correspond to the left-hand and 
right-hand peaks in D3 respectively). For the distribu- 
tion D nsc the exponent r„ sc can be larger than 1 since 
this distribution may extend into the small s region and 
convergence of the integral in (|38|l cannot be ensured. 

Fig. El shows a 3d view of the collapses correspond- 
ing to D\ (sL~ d t ', u 2 L 1 ^ u ^ ■ The lines show three cuts 
of the scaling surface corresponding to u 2 L x l v — 1.21, 
u 2 L 1 / u = and u 2 lM v = -0.56. The collapses of the 
curves corresponding to the different sizes are satisfactory 
within statistical error. The only free parameter in this 
case is df. The best estimation renders a fractal dimen- 
sion df = 2.78 ± 0.05 for such ld-spanning avalanches. 
Similar behaviour is obtained for D 2 (sL~ df , u 2 L x l v } . 
Although, in principle, we have considered df as a free 
parameter, the best collapses are obtained with the same 
value df — 2.78 as that obtained for the ld-spanning 
avalanches. 

The analysis of D 3c and is more difficult. Ac- 
cording to the corresponding distribution relation (see 
Table [HJl, and assuming the scaling hypothesis i|18|) . 120|) 
and one can write: 

N 3 D 3 {s;a,L) = 

= L^'N^uL'/^D^ (sL-d^uL 1 '^ + 
+L- d3 -N i _{uL 1 /")b^ (sL- d3 - ^uL 1 '^ , (39) 

where we have taken into account the fact that for the 
subcritical 3d-spanning avalanches r^d- — 1 and they 
have a fractal dimension rf3_. Although it is possible to 
conceive a DFSS treatment to separate the two contribu- 
tions in H39fl . the hard numerical effort needed as well as 
the associated statistical uncertainties make it very diffi- 
cult. In the next section we will show that it is enough 
to analyze the behaviour of the fc-moments of the distri- 
butions to obtain the critical exponents. 




FIG. 12: Collapses corresponding to f»i [sL~ d f , u 2 L 1/u j . 

The three cuts of the scaling surface are taken at U2L 1 '" = 
1.21, u 2 L 1/v = and u 2 L 1/v = -0.56. Note that on each 
cut we have plotted 5 lines (with different shades of grey) 
corresponding to L = 8, 10, 12, 16 and 24. Symbols on the 
horizontal plane show the behaviour of the first moment of 
the distribution according to the legend in figure |U The line 
is a guide to the eye. 



VII. SCALING OF THE fc-MOMENTS OF THE 
DISTRIBUTIONS 

Besides the scaling of the entire distributions 
D a (s;a, L) that exhibit large statistical errors, it is also 
useful to analyze the behaviour of their fc-moments. For 
those distributions for which the integral in Eq. (|38l) con- 
verges (and, therefore, r a — 1), we can check the corre- 
sponding scaling functions. By using a similar argument 
as that used for deriving equation (|3*5|l . we get: 

L 3 

{s k ) a {a,L) = J2s k D a (s;a,L) = L kd ^ k a (ul}'») . 

s=l 

(40) 

As an example of such collapses, we have indicated the 
behaviour of the scaled first moment of the distribution 
D\ (s; a, L) on the horizontal plane of Fig.^J In this case 
the collapses are obtained without any free parameter. 

As will be seen later, it is more convenient to analyze 
the scaling behaviour of the products N a (s k ) a . By using 
equations (|T3|l and (|10l . one gets: 

N a {a,L){s k ) a {a,L) = L 0+kd «N a {uL 1 '") {uL 1 ^ 

(41) 

Fig.ll3lshows the collapses corresponding to the first mo- 
ment (average size) of D\ and No free parameters 
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are used in this case. Similar scaling plots can be ob- 
tained from the analysis of the second moments with the 
same set of scaling exponents. 
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FIG. 13: Collapses corresponding to Ni(a, L)(s k )i(a, L) (a) 
and N2(a, L)(s k )2(o', L) (b). Symbols indicate the system 
sizes according to the legend. 



FIG. 14: Collapses corresponding to 

Ns-iuL^^^iuiL 1 ^) (a) and N^^L 1 '")^^ 1 '") 
(b). Symbols indicate the values of L\ and L2 used for the 
DFSS analysis according to the legend in Fig. 1101 The inset in 
(a) reveals the power-law behaviour Na-'^\_ ~ (\u2\L 1 / v ) ,33 ~ 
with /3 3 - = 0.024. 



As regards the scaling of N 3 (s} 3 , multiplying expres- 
sion Ij39(l by s, summing over the whole s-range, and 
imposing condition (|37|) . one obtains: 

N 3 (s) 3 = L e+d fN 3c (uL 1 '") *i c (uL l ' v ) + 

+ L d3 -N 3 _ ( uLl/u ) *3- ( MLl/ ") ' ( 42 ) 

This equation can be separated by a DFSS analysis. 
Figs. I14f a) and I14f b1 show the collapses corresponding 
to iV3 C \E , 3 C and N 3 -^>^_ respectively. The only free pa- 
rameter in this scaling plot is the fractal dimension of 
the subcritical 3d-spanning avalanches. The best value 
is d 3 - = 2.98 ± 0.02. Note that the shape of the scal- 
ing function in Fig. I14f bl indicates that, in the thermo- 
dynamic limit, the critical 3d-spanning avalanches only 
contribute to the first-moment for a — a c 



On the other hand, the shape of the scaling function 
in Fig. ll4f a') indicates that, in the thermodynamic limit, 
the subcritical 3d-spanning avalanches may contribute to 
the hrst moment in the whole 112 < range. Note that, 
as revealed by the inset in Fig. I14f a1. the behaviour 
in the region of negative values of U2L 1 / 1 ' is N 3 -^_ ~ 
2 |£ 1/ 'T 3 - with (3 3 - = 0.024 ±0.012. This numerical 
value is compatible with the equation: 



= i>(3 - da_) 



(43) 



This hyperscaling relation, when introduced into equa- 
tion (|42|l . results in a second term that grows with L 3 . 
As will be analyzed in the next section, such a term will 
be responsible for the order parameter behaviour in the 
thermodynamic limit. 
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The analysis of the moments of the non-spanning 
avalanches presents extra difficulties, as occurred in the 
analysis of their number. The expected behaviour is: 

N ns (a,L)(s k ) ns (a,L) = 

+L 3 iV„ s o(a)^ s0 (a) . (44) 

As explained previously, the DFSS cannot be applied, 
given the different dependence on uL 1 / 11 and a of the two 
terms in (|44|) . The possibility of using a trial function is, 
now, more difficult since we cannot make a straighfor- 
ward hypothesis on the shape of ty k sc - In order to fit the 
value of r nsc and d nsc we can analyze the dependence of 
the fc-moment (for k = 2 and k — 3) and its derivatives 
with respect to a at a = a c (u = 0). Data is shown 
in Fig. I15t a) and Fig. I15f with log-log scales. The al- 
most perfect power-law behaviour for different values of 
k and for the derivatives, indicates that the second term 
in (|44Jl plays no role in a c . This is because the expo- 
nent of the first term is much larger than 3. Indeed, 
the best fits are obtained with d rlsc = df = 2.78 ± 0.05 



and t„ sc = 1.65 ± 0.02 which render large values of the 
exponent of the first term (> 5.8). Similar fits can be 
obtained from higher moments with the same values of 
the exponents d nsc and r„ sc . 



VIII. MAGNETIZATION DISCONTINUITY 

In this section we discuss the behaviour of the discon- 
tinuity Am in the magnetization of the hysteresis loop. 
We would like it to behave as an order parameter. For 
large systems, it is clear that only spanning avalanches 
may produce a discontinuity in the magnetization. We 
can evaluate the total average magnetization jump Am s 
due to the contribution of all the spanning avalanches 
(Id, 2d, 3c and 3—). 



Am s = jjN s (s) s . 



(45) 



Fig. llCf a) shows the behaviour of Am s versus a for dif- 
ferent system sizes. According to the scaling analysis in 
the previous section, Am s will behave as: 



Am s (a,L) = 2^L e+d f- 3 [^(ui 1 ^)^! (uL 1 '^ + N 2 (uL 1 ^ * 2 (uL 1 ^ + N 3c # 3c ( uLl/ ") 

+ L d3 — 3 N 3 _ (uL l ' u ^j * 3 _ (uL 1/v )\ • 



(46) 



r 



This equation tells us that Am s will display a mixed scal- 
ing behaviour. The first term in (|46(l accounts for the 
contributions of the ld-spanning, 2d-spanning and criti- 
cal 3d-spanning avalanches. We can define an exponent 
(3 C so that: 



• + d/-3) 



This relation is the same relation that other authors 
have called "violation of hyperscaling" [TTL fT^ . l2^| . From 
our best estimations of 9, v and df, we obtain (3 C = 0.15± 
0.08. 

At this point, it is interesting to compare equations 
()43J) and l|47l ) We would like to note that we could also 
have introduced an exponent 9' that would transform 
Eq. 143fl into an equation similar to (|47|l . Nevertheless, 
the quality of the scalings of the numbers of 3d-spanning 
avalanches in Fig. ^| shows that such an exponent 9' is 
either zero or very small. Moreover, an analysis of the 
behaviour of N%- for uL}l v — ► — oo reveals an exponen- 
tial drift versus = 1 which reinforces the idea that 
there is no need for an hyperscaling exponent 9'. Note 
that a value 9' > implies that the number of subcrit- 
ical 3d-spanning avalanches (3—) will be infinite, in the 



thermodynamic limit. On the other hand, our assump- 
tion that 9' = indicates that behaves as a step 
function in the thermodynamic limit. 

By inserting equations l|43ll and (|47|l into (|46|l one can 
easily read the mixed scaling behaviour of Am s : 



(47) Am s {a,L) L-^'^iu^l^+L-^-^&iuL 1 ^) (48) 



where, (3 c /v = 0.12 and Ps^/v = 0.02. $ is a scaling 
function with a peaked shape (it corresponds to the sum 
of the scaling functions in Figs. I13r a). I13f bl and lMf b)) 
and is twice the scaling function in Fig. I14f b). Con- 
sequently, in the thermodynamic limit, only the second 
term associated to the subcritical 3d-spanning avalanches 
will contribute to the magnetization jump (order param- 
eter). For finite systems, the first term may affect the 
scaling of the data close to a c given the peaked shape of 

This behaviour can be observed in Figs. I16f b') and 
I16f c). where the two possible scalings show the break- 
down of the collapse for U2L 1 / 1 ' < when using the 
exponent (3 c /v and the breakdown of the collapse for 
U2L 1 / 1 ' ~ when using the exponent /?3_/t/. The larger 
the system, the better will be the data collapse in Fig. 
I16r c'l and the worse will be the collapse in Fig. Ilfcif b) . 
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FIG. 15: (a) Behaviour of N ns (a c , L)(s fc ) n s(o" c , L) as a func- 
tion of L for k = 1 and k = 2 in log- log scale, (b) Behaviour of 
the derivative with respect to a of the same two magnitudes. 
In both cases, the lines show the best fits of Eg 1441 and its 
derivative at a = a c . 

IX. DISCUSSION 

Tabic IIIII shows a summary of the exponents that 
characterize the avalanche numbers and distributions ob- 
tained from our numerical simulations. We would like to 
point out that such exponents are independent of a and 
L in a very large region around the critical point both for 
a > <t c and a < a c simultaneously. Such an achievement 
has not been possible in previous analyses, even with 
larger system sizes. The reason is that some of the con- 
tributions we have identified (namely 3— and nsO), which 
reduce finite size effects, were previously neglected. 

In Tablc lTTTl we also indicate previous estimations of the 
exponents found in the literature |13|. The comparison 
is quite satisfactory. Let us analyze the eight exponents: 

• Although the value of v does not fall within the 
error bars in Ref. Il3l we have already argued that 
the exact definition of the scaling variable u used 



exponent 


best value 


values in Ref. 13 


V 


1.2 ±0.1 


1.41 ±0.17 


e 


0.10 ±0.02 


0.15 ±0.15 




2.02 ±0.04 




df 


2.78 ±0.05 


2.98 ±0.43 (= 1/av) 


d 3 - 


2.98 ±0.02 






1.65 ±0.02 




0c 


0.15 ±0.08 




/3 3 - 


0.024 ±0.012 


0.035 ± 0.032 (= (3) 



TABLE III: Summary of the values of the exponents obtained 
from the simulations in this work. We have indicated the 
names of the exponents from Ref. llH . whose definition does 
not exactly correspond to our nomenclature between paren- 
theses. 



for the collapses may introduce some deviations in 
this value. By using u\ we obtain v — 1.14 and 
using U3 we obtain v = 1.4. 

• As regards 9 our value is in agreement with the 
value previously reported [13j (We would like to 
note that in Ref. Il3l the authors also report a 
value of0.015±0.015 probably due to a misprint). 
The fact that this exponent is non-zero implies that 
there are infinite spanning avalanches at the critical 
point in the thermodynamic limit. 

• As regards 9 nsc , to our knowledge there are no pre- 
vious finite-size scaling analyses of the number of 
non-spanning avalanches. 

• Concerning df and d%-, the numerical values are 
consistent with the value df — 2.98 ± 0.43 esti- 
mated previously . We shall note that this pre- 
vious estimation was obtained from the analysis of 
the distributions of non-spanning avalanches. It 
should therefore correspond to our df and not to 
our c?3_ (which corresponds to the subcritical 3d- 
spanning avalanches). Note also that the difference 
between df and dz- suggests that there might be 
real physical differences between such two kinds of 
avalanches. The possibility of distinguishing them 
in the numerical simulation will be studied in a fu- 
ture work. 

• The exponent r nsc , according to our definitions, de- 
scribes the scaling behaviour of the distribution of 
critical non-spanning avalanches. Previous mea- 
surements of a similar exponent, have analyzed 
N ns without distinguishing between critical (nsc) 
and non-critical (nsO) non-spanning avalanches and 
have not considered the fact that the system is fi- 
nite. We can estimate what the value of an effective 
exponent r e ff will be for the distribution of non- 
spanning avalanches for very large systems. From 
equation 144|) . taking k > 1, and large values of L, 
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only the first term in the sum survives, so that: 
N ns (a)(s k ) ns (a) = 



_ j j e nsc +{i+k- 



(49) 

On the other hand, in the same limit, the analysis 
of equation l|26|) renders: 



N ns (a) =L 3 N ns0 (a) 



(50) 



Combining the last two equations, we get an es- 
timation for the pseudo-scaling behaviour of the 
fc-moment of the non-spanning avalanches: 



= L 



(51) 



If one approximates N ns o(cr) by N ns o(a c ) and im- 
poses (s k ) ns ~ L-^ff- k -^ d fS k (uL x l v ) it is pos- 
sible to deduce that the effective exponent is T e f t = 
T nsc + (3 — Onsc) /df. From our numerical estima- 
tions of the different exponents in Table IIIII one 
obtains T e ff — 2.00 ± 0.06. This value is in very 
good agreement with the value r e ff = 2.03 ± 0.03 
found previously Nevertheless, we would like 
to point out that according to our analysis, such an 
exponent is not a real critical exponent and, there- 
fore, will depend on a for finite systems as has been 
found previously |16| . 

• As regards the values of (3 C and /?3_ we would 
like to note that previous analyses have not iden- 
tified the two contributions to Am s . It is there- 
fore not strange that different values have been ob- 
tained previously: 0.17 ± 0.07 @, 0.0 ± 0.43 [HI, 
0.035 ±0.028 |ig. The larger the system, the closer 
the effective exponent becomes to /?3_. 

Finally, it is interesting to compare the behaviour of 
spanning avalanches, with the problem of percolation 
|23|. In percolation, the number of percolating clusters 
behaves as a step function, in the thermodynamic limit 
for d < 6, exactly as Ns_. The order parameter is, in this 
case, the probability for a site to belong to the percolating 
cluster. However, this is precisely what we are evaluating 
by the function jV3_(s)3_/L 3 which is the second term 
in (|46|l and is the only relevant term in the thermody- 
namic limit. As occurs in percolation, the hyperscaling 
relation Q43JI between (3^-, v and d-$- is fulfilled since 
only one infinite avalanche contributes to the order pa- 
rameter for er — > a c from below. Moreover, in the perco- 
lation problem for d > 6 the number of percolating 
clusters exhibits, besides the step function, an extra 8- 
function singularity at the percolation threshold. In our 
case (3d-GRFIM) we also have such a contribution at 



a = (T c , which we have identified as Id, 2d and critical 
3d-spanning avalanches [the first term in (|46|) ]. The exis- 
tence of such an infinite number of avalanches exactly at 
a c (the number of which grows as L e ) implies the break- 
down of the hyperscaling relation j3 c = v[S — (9 + df)]. 



X. SUMMARY AND CONCLUSIONS 

In this paper we have presented finite-size scaling anal- 
ysis of the avalanche numbers and avalanche distribu- 
tions in the 3d-GRFIM with metastable dynamics. After 
proposing a number of plausible scaling hypotheses, we 
have confirmed them by obtaining very good collapses of 
the numerical data corresponding to systems with sizes 
up to L — 48. 

The first result is that, in order to obtain a good de- 
scription of the numerical data, one needs to distinguish 
between different kinds of avalanches which behave dif- 
ferently when the system size is increased. Avalanches 
are classified as being: non-spanning, ld-spanning, 2d- 
spanning or 3d-spanning. Furthermore, we have shown 
that the 3d-spanning avalanches must be separated into 
two classes: subcritical 3d-spanning avalanches with 
fractal dimension e?3_ = 2.98 and critical 3d-spanning 
avalanches with fractal dimension df = 2.78, as the Id— 
and 2d— spanning avalanches. Non-spanning avalanches 
occur for the whole range of a. We have also proposed a 
separation between critical non-spanning avalanches and 
non-critical non-spanning avalanches in order to obtain 
good finite-size scaling collapses. The non-critical non- 
spanning avalanches arc those whose size is independent 
of the system size and whose number scales as I? . The 
critical non-spanning avalanches also have a fractal di- 
mension df — 2.78. 

The second important result, is the scenario for the 
behaviour in the thermodynamic limit: below the crit- 
ical point, there is only one subcritical 3d-spanning 
avalanche, which is responsible for the discontinuity of 
the hysteresis loop. Furthermore, at the critical point 
there are an infinite number of Id-, 2d-, and 3d- critical 
spanning avalanches. 

For finite systems, the six different kinds of avalanches 
can exist above, exactly at and below a c . The finite- 
size scaling analysis we have performed has also en- 
abled us to compare different scaling variables u, which 
measure the distance between the amount of disorder 
in the system a and the critical amount of disorder 
a c . The best collapses are obtained using the variable 
ii2 = (<t — <t c )/(T c + A [(a — <r c )/cr c ] 2 with A = —0.2. 

So far the analysis presented in this paper, is re- 
stricted to the analysis of the numbers and distributions 
of avalanches integrated along half a hysteresis loop. Our 
analysis of the average magnetization discontinuity Am 
starts from the hypothesis that only spanning avalanches 
may contribute to such a discontinuity. However, as a 
future study, we suggest that the measurement of cor- 
relations in the sequence of avalanches and the analysis 
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of non-integrated distributions, may reveal details of the 
singular behaviour at the critical field H c . For instance, 
non-spanning avalanches could show a tendency to accu- 
mulate in H c , in the thermodynamic limit. This could 
change some of the conclusions reached in this work. 

As a final general conclusion we have shown that it 
is not necessary to simulate very large system sizes to 
estimate the critical exponents for this model. In order 
to identify the different kinds of avalanches it may even 
be better to analyze small systems with larger statistics. 
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FIG. 16: (a) Behaviour of the total contribution of the span- 
ning avalanches to the magnetization jump as a function of 
a. The inset reveals the crossing behaviour at a < a c . (b) 
Scaling of Am s by considering the Id, 2d, and the critical 3d- 
spanning avalanches. Note the lack of collapse for the region 
uiL x l u < 0. (c) Scaling of Am s by considering the subcritical 
3d-spanning avalanches. Note the lack of collapse for the re- 
gion ~ 0. Symbols indicate the system sizes according 
to the legend. 



